A nonlinear hydrodynamical approach to granular materials 
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We propose a nonlinear hydrodynamical model of granular materials. We show how this model 
describes the formation of a sand pile from a homogeneous distribution of material under gravity, 
and then discuss a simulation of a rotating sandpile which shows, in qualitative agreement with 
experiment, a static and dynamic angle of repose. 



I. INTRODUCTION 



The nature of the theory which describes the macroscopic transport in granular materials remains an unsolved 
problem. We propose here a candidate theory, inspired by features of nonlinear hydrodynamics (NH). Given its success 
in treating transport in a variety of complex systems NH seems natural to use in this context. In developing any 
continuum model of granular materials, however, one is challenged by the need to incorporate some remnant of the 
(4!^ ' discrete nature of the underlying material. 
O [ The development of our theory is based on the hypothesis that the states of a system can be specified in terms of 
a few collective variables (at sufficiently long length scales), and that these states are connected in time by the local 
' conservation laws supplemented with constituitive relations. Typically, the collective variables used are the conserved 
densities which, in the simplest view of granular materials, are conservation of mass and momentum. Conservation of 
energy is more complicated in this case than in simple fluids, due to locally inelastic processes which may transfer 
'■^ , energy into internal degrees of freedom. Although we know how to include the energy density into the description, 
^ ■ we begin with a more primitive theory which uses only the mass density p and the momentum density g. 
O I The key difference between simple fluids and granular materials is that fluids organize over short time scales 
to be spatially homogeneous, while granular materials can exist over very long times in spatially inhomogeneous 
metastable states: metastable, because individual grains do not always pack together efficiently; long-lasting, because 
1: the large masses of the grains prevent thermal fluctuations from adjusting particles into a tighter configuration. Any 
^ ■ complete description of granular materials must be able to explain the nature of these quenched inhomogeneous states. 
Experiment suggests that a sandpile contains "force chains" |^ which serve to support the pile in the presence of 
external stress; these chains are surrounded by regions of sand which are comparatively unstressed, and which allow 
] the pile to flow and to settle under vibrations. 
■ Since these sandpiles (in the absence of forcing) are static, the momentum is negligible, and we are left with the 
' density field as our sole tool in describing the pile. Our primary hypothesis is that the density field alone is sufficient 
. to capture the metastable nature of the system when at rest. Specifically, we describe the force chains as regions of 
sand in which the packing of the grains is ideal, having the maximum possible density. We will call these regions 
"close-packed" . The other, "loose-packed" regions have their grains arranged in some non-optimal configuration, with 
y • a slightly smaller density. We propose that these small variations in the density field are sufficient to describe the 
r metastable nature of sand. 

(— I Using dynamical equations based on local conservation of mass and momentum, this quenched pattern in the density 

Q can be prepared in a dynamic process driven by an effective free energy with minima characterizing the loose-packed 
O and close-packed regions. This presents the opportunity to grow a sand pile under the influence of gravity starting 
^ from an initial homogeneous state, and in the construction of the sand pile one builds up, in a self-consistent manner, 
• I— I . an inhomogeneous equation of state characterized by a nontrivial stress tensor. 

' The basic structure of our theory, based on local conservation laws, is straightforward. More difflcult is the inclusion 
^ ] of these competing close and loose-packed states. Clearly these require a detailed description of the shortest length 
- - ' scales treated in the model; there is a competition of length and energy scales of a type not encountered in simple 
fluids. Such a short-range description will depend on details of the particular granular material; we will be satisfled 
with a theory that demonstrates some generic properties of sand, as described in the next section. It seems reasonable 
that one may then be able to refine the model to include variations in the shapes and types of individual grains. 

While the short-range length scales give difficulty, the large-scale nature of our model presents us with the advantage 
of a description close in scale to those phenomena which are most visible to experiment and casual observation. One 
ultimate goal is to address the existing macroscopic shaking and rotation experiments. Also, such a theory should 
eventually be useful in analyzing the surface states produced in large-amplitude shaking experiments 

There have been a number of attempts Q to use hydrodynamics to understand the dynamics of sand piles. These 
efforts differ from the one advanced here in that they do not propose to follow the full evolution of the system. There 
has also been a significant theoretical effort to describe granular materials through the use of kinetic theory p3| . This 
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approach is organized at a more microscopic level, where the connection to the macroscopic, static, metastable phase 
of the system is less direct ||l^ and coupling to experiment more difficult. In many of these cases it is necessary to 
"liquefy" the granular material in question, giving it so much energy that it loses its metastable quality which seems 
to be crucial in describing certain aspects of a sandpile's behavior. 

In this paper we begin by describing some generic features of sand which we hope to incorporate into our model. 
We then proceed to a development of our dynamical equations from a generalized Langevin equation. This will lead 
into a discussion of the choice of a driving free energy, which will be crucial in creating the quenched pattern in the 
density field. 

To test the validity of the theory we turn to simulation on a two-dimensional lattice. We begin with a rectangular 
box containing a homogeneous distribution of sand under the influence of gravity, and we show that it does indeed 
form a sand pile with an inhomogeneous distribution of loose and close-packed states. We then turn to simulations in 
a circular container, first forming a sand pile using the same technique, and then rotating the pile at different speeds. 
We compare the resulting angles of repose and oscillations about those angles to corresponding experimental data in 
the literature. 



II. PROPERTIES OF GRANULAR MATERIALS 



We want to construct our model to be compatible with the following aspects of granular materials: 

1. Granular materials have a clumping property; simulations of simple systems of inelastically colliding balls provide 
evidence for this phenomenon. One possible explanation for this clumping is the theory of inelastic collapse 

when two particles collide inelastically, they lose energy from their translational degrees of freedom and 
so recede more slowly than they approach each other. On average, the particles stay closer together than if 
they had collided elastically, and regions of higher density build up. This can also be described in terms of a 
hydrodynamical instability [|l5| . 

2. Thermal energies in a sand pile are very small compared to, for instance, the average gravitational energy. Thus 
we should be able to ignore thermal noise. 

3. Because of the large masses of the grains, one does not expect a significant vapor pressure above the interface 
as is found in liquid-gas systems. Indeed we expect a very dilute gas of grains above the pile, whose density 
exhibits a Boltzmann distribution due to the effect of gravity. 

4. Sand piles are strongly driven by gravity and, because of the lack of thermal effects, the dense sand pile is 
separated from the dilute gas above by a sharp interface. 

5. Granular materials are strongly disordered, with metastable structures forming upon creation of a pile. This 
is also due to the lack of thermal effects. As mentioned in the introduction, there is evidence of stress chains 
running through the bulk of sand piles, which may be associated with the ramified clumps formed in the more 
dilute systems studied using kinetic theory. The important point for us is that there is competition between 
somewhat more dense domains (stress chains and arches) and other, more loosely-packed regions in the pile. 

6. Granular materials are stiff: when moved or jostled the pile can for a time behave as a solid object. Thus the 
sand pile should have a relatively uniform density. 

7. One specific example of this stiffness is that a sand pile can maintain a non-horizontal surface. For example, 
when a pile of sand in a container is rotated about a horizontal axis, it does not flow immediately but waits 
until its surface passes a certain critical angle (called the angle of repose) with the horizontal. According to 
Reynolds p6) , this effect is due to the static interlocking grains that make up a granular material: the pile must 
dilate first, giving these grains freedom to move about, before flow is possible. 

8. It is well known that granular materials can be very sensitive to boundary and finite-size effects; the famous 
Brazil-nut-to-the-top phenomenon [pT| , for instance, is strongly infiuenced by both. 



We will touch on these points while developing our model. 
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III. DYNAMICAL EQUATIONS 



Inspired by nonlinear hydrodynamics, we begin with the generahzed Langevin equation |18, [L9| 

(Throughout this section, summation is imphed over aU repeated indices.) The variables ipa are the "slow" fields of 
interest in the problem: as mentioned earlier, these are the density field p{x) and the momentum field g{x). The first 
term Vq.[?/'] is the streaming velocity, and corresponds to the reversible terms in typical hydrodynamical equations. As 
we will see, it depends on the Poisson brackets of the slow fields, as well as the derivatives of the Landau-Ginzburg- 
Wilson (LGW) effective free energy F. The second term is dissipative in nature, and TafS is the symmetric matrix of 
dissipative coefficients. Because granular materials are essentialUy zero-temperature systems (property there is no 
thermal noise driving the system. 

For our system, the free energy F can be broken up into kinetic, potential, and external parts: 

F = Fk+Fp + Fe. (2) 

From thermodynamics, one has the result that the variable conjugate to the momentum density is the velocity field: 

SF _ _ gijx) 

'"i[x) = —F^i (3) 



Sgi{x) p{x) 

where i is a vector label. The kinetic energy contribution to the free energy thus has all of the momentum dependence: 

while Fp and Fe are functions only of the density. We will make further assumptions as to the form of Fp below. 
Fe is the free energy due to external forces: for instance, in the presence of a uniform gravitational field 

Fe= I d^xgp{x)iz-zo), (5) 



where the scalar g is the acceleration due to gravity and zq is the bottom of a confining box. 
The streaming velocity Va (in Eq. (^) is given by the equation 

SF 

V4^] = {^^,i:0}— (6) 

where the indices a, f3 run over the set {p,g}- We calculate the Poisson brackets by identifying the fields ipa with 
microscopic variables, evaluating their Poisson brackets, and expressing the results in terms of the ijja, getting the 
standard results 

{p{x),g,{x')} = ^V^[S{x-x')p{x)] (7) 

and 

{g,{x),g,{x')} = -Wj[d{x^x')g,{x)] + W'Mx - x')g,{x)]. (8) 
With these and the appropriate derivatives of the free energy, we can calculate the streaming velocities: 

Vp = -V-g (9) 

and 

By choosing Tpp = 0, one easily obtains the usual continuity equation 
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The components of the damping tensor associated with the momentum density, 

^'tj = Tgigj, (12) 

can be expressed in terms of the viscosity. If L were independent of the fluctuating fields, then we could write it in 
the general form 

Lij = ~Vii,kj'^i^k (13) 

where 

is the viscous tensor, and 770 and (^q are the "bare" viscosities. It is more realistic, however, for the dissipation to 
depend on the density of the sand, and so we introduce a function 4i{p) into the damping tensor: 

Lij = -rjiiM^i {<i>{p)^k) ■ (15) 

We shall say more about (f>{p) below. It should be emphasized that, in a system with strong nonlinearities and spatial 
inhomogeneity such as ours, there exists no simple relationship between these bare viscosities and their physical 
counterparts. 

We can now write the remaining Langevin equation: 

at ~ '-''Sgj 

= -pVj-^ - pgS^^ - Vjig,gj/p) - (16) 

At this point, we must choose a form for Fp. In this simple theory, we make the assumption that Fp is of the 
square-gradient form 



Fp= 



fip) + IciVpf 

where f{p), the free energy density, is a local function of p, and c is a positive constant pO|. Thus we have 



(17) 



-P^t-Q^ + cpWiV'^p - VJ{g^gJ/p) + VijM^j [HpWk{gi/p)] - gpSiz- (18) 

It should be noted that we can rewrite this equation as the divergence of a stress tensor plus a term representing 
external forces: 



^ = -V,a,+Ff. (19) 



This is the usual continuity equation for conservation of momentum, just as our expression for the density's evolution 
is the continuity equation for conservation of mass. 

For calculational reasons, we prefer to work not with the momenta g, but the velocity fields v, defined by Eq. (|^) 
above. In terms of the density and velocity fields, the equations of motion become 

^ = -y-{pv) (20) 

and 

dvi df 1 

-g^ = ^'^^■^ + cViV^p - + -miM^j [4'{p)^kVi\ - gSiz. (21) 

To proceed further, one needs to choose forms for (j){p) and the viscosities. For very low densities we expect that the 
dissipative part of the equation should vanish; for this to happen, (j>{p) must go to zero at least as fast as p^. Higher 
powers of p seem to make our numerical calculations more stable, so we have used ip — P^- This issue is not very 
important in the dense sand pile where p « 1. We choose Co ~ f ^0 in Eq. ( p^ ) so that rjij ^i is isotropic: 

VijM = VoiSijSki + SikSji + SiiSjk) (22) 

In this theory, the input most crucial in specifying sand-like behavior is the free energy density /(p). With an 
appropriate choice of / we will be able to create a system that is very different from conventional fluids. 
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IV. CONSTRUCTING THE FREE ENERGY DENSITY 

We will construct the free energy f{p) using three steps, as outlined in Figure The mathematical details of the 
description of the free energy can be found in the Appendix. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 

P P P 



FIG. 1: Three steps to our free energy density: a) with no met ast ability; b) with a barrier to create metastability, and c) with 
wells to stiffen the material. 

We can think of /(p) as a potential where the system picks out values of the density which correspond to minima 
of /. A key assumption is that, from Property |l|, grains of sand tend to clump together, and so the potential has a 
minimum which we arbitrarily place at p = 1. The grains themselves are incompressible, so the potential rises rapidly 
for densities higher than p = 1. Similarly, the potential becomes large as one approaches zero density, to prevent 
unrealistic negative values of the density anywhere in the pile. 

Calculations show that this potential, so far, is enough to capture some properties of sand: the sharp interface of a 
sand pile, for instance. However, it is still missing one key element, and that is the metastability found in experiment. 
Any pile formed using this free energy would be homogeneous in density, and show no evidence of the close-packed 
and loose-packed regions mentioned in property |^. Without these regions, our material would be little more than a 
slightly compressible liquid. 

In Figure |^b, we introduce this metastability into the system by placing a barrier into the potential just to the 
left of its minimum. We now have the following picture: the minimum represents the optimal packing for the sand 
grains. As sand comes together, it increases the density of the pile until it reaches the density of the barrier where 
it is frustrated: it can no longer proceed to the optimal packing. However, as the sand begins to pile up, pressures 
in the pile can provide enough of a push so that some regions may pass the barrier and reach the minimum of the 
potential. 

This potential has all of the right features, but in simulations it proves to make a sandpile that is too soft and 
liquid-like. To stiffen the pile, we introduce two wells, as in Figure ^: one for the close-packed states and one for the 
loose-packed. With the wells, this is the free energy we have used in our calculations. Future possible modifications 
to this free energy will be discussed in our conclusion. 

V. NUMERICAL ANALYSIS 

Because the proposed model has a strongly nonlinear nature, we are forced to study this system's evolution numer- 
ically, hoping that this will inspire subsequent analytical work. The first order of business is to determine numerically 
whether this system of equations, running forward in time, can prepare a state that looks like a sand pile. In principle, 
this should be straightforward to do. In practice, there are a number of difficulties in putting these partial differential 
equations onto a space-time lattice. 

1. Lattice spacing: It is understood that, to solve PDE's numerically, one must keep the time step much smaller 
than the lattice size; otherwise, the system develops physical instabilities. We have chosen a typical convention 
of setting Ax = 1 and At = 0.001. 

2. We know that there are important finite-size effects in granular systems, and so we must be prepared to 
examine our model in various sizes and shapes of containers. This paper focuses on two such containers: a 
rectangle 200 units high and 100 units wide, and a circle with a 100 unit diameter. Other containers have been 
tried as well, yielding similar results. 



6 



3. One of the discouraging aspects of this approach is the large number of parameters governing the system. There 
are parameters to control the intrinsic free energy, the square-gradient energy, as well as those for external forces 
like gravity. There is the finite size of the box and the nature of the boundary conditions. Dynamically, there 
are the viscosities and the initial level of the density and velocity fields. In this paper, the values of these 
parameters have been chosen to obtain structures that resemble sand piles. Most of these choices are reported 
in the appendix where the form of the free energy is described in detail; in addition, we take c — 10 and 770 — 12. 
Further work is needed to map out the range of parameters for which we obtain physical behavior, and determine 
which parameters are ultimately important. 

4. This numerical problem would be rather straightforward if not for a set of persistent unphysical instabilities; 
without care the system can and will explode in a very unnatural fashion. These explosions are of two, closely 
related types: runaway velocities and negative densities. To alleviate the first, we have taken the pragmatic 
step of locally averaging the velocity if it gets larger than some empirically determined cutoff. 

Our free energy density has been designed with a barrier at p = 0, but negative densities do tend to show up in 
our calculations, specifically in the dilute region above the sand pile. We take two steps to counter them: when 
the density of a site is small and positive {p < 0.05), we average the site's velocity with half of the neighbors' 
average velocity. In addition, when the density dips below zero, we bring in density from its four neighbors to 
bring it above zero. In treating these instabilities, we have been very careful not to violate conservation of mass. 



VI. FORMING A SAND PILE 



As a first example of how our model works, we begin with a nearly uniform distribution of sand which, under 
the influence of gravity, forms a sand pile at the bottom of the container. Our container is a two-dimensional box, 
overlaid with a square lattice 100 units wide and 200 units high. Boundary conditions on the walls of the container are 
non-slip. Our initial conditions are p = 0.5 + Ap and v — Av, where Ap and Aw are both small (0.001 rms) random 
perturbations in the system. Because mass is conserved and because the density inside the formed pile should have a 
density about 1, we expect the interface to form in the middle of the box, around z — 100. In principle, one should 
average over ensembles of initial conditions. That is not necessary for our purposes here, but will become important 
in more quantitative work. 

The first question is whether the model captures the overall dynamics of the system: our physical picture of the 
situation has the sand particles all under the influence of gravity, so at first there will be a net acceleration. As 
particles begin to hit the bottom of the container, their velocity drops to zero, and the system as a whole decelerates 
until it ends up at rest. To see if this is true we measure the average kinetic energy 

where the sum is over the entire box, and V is the total number of lattice sites. Figure |^ shows a plot of the kinetic 
energy for two runs with different gravitational accelerations. 

In both cases the behavior is just as expected, with an immediate acceleration followed by a deceleration to rest. 
Also interesting to note is the difference between the curves: the run with stronger gravity causes the pile to form 
faster, with a higher maximum kinetic energy resulting from the higher initial potential energy of the system. 

To get a clearer picture of the stationary state our system is settling into, we consider the density profile: 

where is the number of sites across the lattice. 

Figure ^ superimposes snapshots of this profile at seven instances in time. Initially, the profile is a flat line 
P(z) = 0.5, representing the initial homogeneous distribution. As time progresses, the density begins to increase 
for low z, indicating pile formation; at the same time, the density at higher altitudes is decreasing. By i = 30 the 
system has settled into a phase-separated system of high-density pile below and low-density "gas" above, with a sharp 
interface between the two (agreeing with property ||in section II.) 

The height of the pile can be defined as that value of z for which the profile reaches P{z) — 0.5, or half the average 
density of the pile [gl] . This is plotted as a function of time in Figure ^, along with the vertical coordinate of the 
system's center of mass: 

= y p(z x) ^^^^ 
2^x,z py^' -^j 
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FIG. 2: Kinetic energies for g = 1 and g — 0.5 
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FIG. 3: The profile of the pile for g — 0.5, t = 5, 10, 15, 20, 30, 40, 50. The variable z represents the height above the bottom of 
the container, which is at z = 0. 



with the sums over all the points in the lattice. 

At first the center of mass moves downward in time parabolically, as it should with the system in free fall. As the 
sand starts to pile up, less and less of the sand is in motion and the center of mass approaches a constant height of 
52.6. The height of the interface moves contrary to this, of course, settling in at 102, about halfway up the container. 
Since the average density of the dense pile should be about 1 , and we began with a uniform distribution of density 
0.5, the position of the interface is just a little bit higher than our expectation that the container be half full of sand, 
suggesting that the pile has an average density slightly less than one. The center of mass is slightly higher than the 
midpoint between the base and the interface of the pile, suggesting that the pile is somewhat top-heavy, which is 
surprising. 

Since on the large distance scale our system looks like a sand pile, we turn our attention to the small-scale features 
inside the pile. Specifically, we want to investigate the competing regions of loose-packed and close-packed states, 
mentioned in property ^ of section II, corresponding to the two wells in the free energy. Recall that these wells are 
separated by a barrier, which we have placed at p = 0.99; for definiteness, we say that all sites with p > 0.99 are 
close-packed and all sites with 0.85 < p < 0.99 are loose-packed. 

In Figure |5| we count the number of loose and close-packed sites as a function of time, for two different strengths 
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FIG. 4: The location of the interface, and the vertical component of the center of mass, for g = 0.5. The interface does not 
really exist until a pile begins to form, so we did not begin measuring the location of the interface until some time after t = 5. 
The sharp peak in the interface's curve occurs when the small bump to the right in figure |^ hits the forming pile below. 
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FIG. 5: Number of loose-packed and close-packed states for g = 0.5 (left) and g — 1. 



of gravity. Naturally, the loose sites start to form first; some of these are then pushed across the barrier and become 
close-packed due to the weight of sand on top of them. Both graphs show the numbers of loose and close-packed sites 
approaching a constant as the pile settles. Notice that the final ratio of low-density to high-density sites is sensitive 
to parameters such as gravity: a stronger gravitational field has a greater capacity to compact the sand and create 
more close-packed sand. 

Figure ^ shows how these loose and close-packed regions are distributed through the pile for the same two values 
of g. The top halves of each pile share similar features: both have a layer of loosely-packed sand on top, where there 
is no other sand to weigh it down and compress it. Immediately below the surface in the center there is a region of 
tightly-packed sand, which dips farther below the surface in the center of the pile than to either side. This region is 
flanked by columns of loosely-packed sand, which are in turn flanked by dense sand up against the edges of the box. 
The differences in the two piles are mostly at the bottom: the run with less gravity is mostly loose-packed in the bulk 
of the pile, while the run with more gravity is tighter. Note however that in neither case is the density uniform: one 
can find patches of both types in both piles. In particular, the patterns in the center of the g = 0.5 pile may hint at 
the stress chains that are seen in experiment. 
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g=0.5 g=1.0 

FIG. 6: Pattern of loose (grey) and close (black) states for g = 0.5 (left) and g — 1. 

VII. ROTATING A SAND PILE 
A. Set-up 

One of the standard ways to probe a granular system is to rotate it about some horizontal axis. This method 
demonstrates several characteristics of a granular material, and the most important of these is the angle of repose. 
Actually, there are (at least) two such angles associated with rotation. If one begins with sand in a cylindrical 
container, having a horizontal free surface, and then begins to rotate the cylinder about its axis, the surface of the 
sand will increase in slope along with the container's rotation until it reaches some critical slope, at which point an 
avalanche will restore the surface to a more horizontal state. This we might call a "static" angle of repose, 9s- If one 
continues to rotate the cylinder, then one of two things will happen: 

1. If the rotation is slow enough, then the flow along the surface will come to a stop before a further increase in 
slope should cause an avalanche again. The result is a periodic, intermittent flow, and the interface oscillates 
about some average angle. 

2. If the rotation is faster, the flow from the first avalanche does not have time to stop before the next avalanche 
begins. The system settles into a state where there is a continuous flow along the surface, and the interface 
maintains a constant angle with the horizontal. 

In both cases, we may take the mean angle that the surface makes with the horizontal over time and call it a "dynamic" 
angle of repose, 9d- In the periodic case, a quantity as important as the angle of repose is the size of the fluctuation 
about that angle, 6d- The actual values of these angles seems to depend experimentally on a number of factors, 
including particle size and shape, the humidity of the air, how the pile is formed, boundary conditions, and so forth. 
The experimental results of Jaeger et al. for instance, flnd that spherical glass beads with diameters of about half 
a millimeter show a dynamic angle of repose of 26°, with a fluctuation of 2.6°; while rough aluminum- oxide particles 
with the same diameter show a higher angle, 39°, with a fluctuation of 5°. These seem to be typical values. What 
will be most important in our qualitative analysis is that these angles exist and are non- negligible. 

We set out to find evidence of these two angles in our model. To better match experiments in this field, we 
move from a rectangular box to a circular one having a diameter of 100 lattice sites and the same no-slip boundary 
conditions. We grow a sand pile here as we did before: Figure ^ shows the pattern of loose and close-packed states in 
a pile formed in our circular container. 

To best mimic rotation under experimental conditions, one would ideally set up rotating boundary conditions, 
giving a constant velocity to the sand at the edge of the container. This has turned out to be difficult to do in 
practice: one reason is that, due to the lattice, our boundary is not a perfect circle, and there is a tendency for mass 
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FIG. 7: Pattern of loose (grey) and close (black) states for </ = 1 in a circular container. It is interesting that this looks 
more like the lower gravity distribution in the rectangular box, in Figure ^ with a clump of high density floating on top of a 
low-density sea. This is probably due to the smaller amount of total mass in this system. 

to "leak" in or out of the boundaries when the boundary conditions are not strictly no-slip. Because of this, we choose 
to implement rotation by rotating gravity with a constant period of rotation T . 

To look for an angle of repose we need to measure the angle that the interface makes with the horizontal (that is, 
the normal to gravity). The most direct way to do this is to fit the interface of the pile with a line, and measure 
the angle that this line makes with the horizontal. This is indeed one of our probes: we define as our interface those 
points which are in the pile (p > 0.5) which have a nearest neighbor outside of the pile, and, to reduce boundary 
effects, we throw out those points that are not within one half radius from the center. 

This method should be ideal, but because of the finite size of our system, this measure ends up depending on only 
a few points which makes it sensitive to various local perturbations. For this reason, we introduce another measure 
which we call the "bulk angle" (where the first might be called the "surface angle"). This latter measure is simply 
the angle that a line passing through the center of the container and the center of mass makes with the vertical. This 
measure depends on the entire system and so is less sensitive to noise. It will be equal to the surface angle in the 
event that the pile has reflection symmetry across the vertical axis; it will differ from the surface angle by some fixed 
amount if the shape of the pile during rotation is asymmetric but constant. Figure ^ shows a plot of the two measures 
for a typical run. 



surface ^ngle 




FIG. 8: The behavior of the surface and bulk angles for a typical run (T = 200). Notice the overall similarity in the shape 
of the two plots, although the bulk angle has a much smoother curve. The bulk angle is consistently above the surface angle, 
suggesting that the pile has an asymmetric shape. (All angles are measured in degrees.) 
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B. Static Angle of Repose 



With the surface angle and bulk angle measures in place, let us look systematically at our data. First, we consider 
the behavior of the pile in the first moments of rotation, with three different periods T (Figure 




angle rotated angle rotated 



FIG. 9: The behavior of the bulk (left) and surface angles during the first 45° of rotation, for three rotation speeds: T — 200, 
280, and 400. The horizontal axis shows the angle through which the container has been rotated. 

In all three cases the behavior is similar: as the container turns, the interface's slope increases, until its angle 
reaches some maximum and begins to decline due to one or more avalanches. Clearly, there is a static angle of repose 
here, which depends on the rotation speed. Notice, however, that if there were no activity in the sandpile before the 
first maximum in the angle, then the curve would follow the straight line denoted in the figure as "solid" until turning 
downward. Instead, it seems that the surface is losing slope, lagging behind the container, from a very early time: 
there is some preliminary flow in the pile, even before the first major avalanche. 

Figure |l^ contains six snapshots of the flow during this initial period, with T = 200. All of the pictures here are in 
the lab frame (that is, gravity's frame), so the flrst picture {t = 4) shows the pile moving in unison with the container, 
as if it were a solid and fastened to the walls. At t = 8, however, things are beginning to change: the center of rotation 
seems to have moved to the right and downward, and the sand coming up on the right side is beginning to curve over 
to the left. The result of this action is seen in the next picture, where there is a definite flow along the surface of 
the pile. Notice that, even though there is surface flow, it is not enough to prevent the angle of the interface from 
climbing, according to Figure ^. The next picture, at t = 16, shows the surface flow continuing, and also that sand is 
beginning to climb up the side of the right wall, due to our no-slip boundary conditions. At t ~ 20, a major avalanche 
begins, so that the angle of the interface begins to fall (see Figure ^; this corresponds to 36° on that figure's horizontal 
axis). The last picture {t — 30) shows a much larger flow of material than seen earlier. 

In short, there is clearly an initial period where the sand moves with the container and where there is no surface 
flow, evidence for a nonzero static angle of repose. Recall from section II that this effect is due to the need of the 
pile to dilate before it can flow. In our model, the close-packed sites are the ones which are restricted in their ability 
to dilate (because of the barrier in the free energy). Thus, if the pile is dilating during this flow, we expect that the 
number of close-packed sites in the pile must be decreasing. Figure |ll| shows that this is indeed so. 



C. Dynamic Angle of Repose 

We next consider how the pile behaves under further rotation with three or more turns of the system. Figure 
above shows the behavior of the bulk and surface angles over three complete revolutions, for T = 200. Figure |l^ 
below compares the bulk angles between T = 200 and T = 1600 over three revolutions. 

Most notable in Figure |l2| is that the average slope of the interface is larger for a faster rate of rotation. Also, the 
slower system deviates much less from a constant angle than does the faster one, with more jaggedness suggesting 
frequent small avalanches. These points suggest an inertial effect: in the faster case the sand does not have as much 
time to avalanche and level its interface, while the slower case has a number of short avalanches that help to keep its 
interface closer to the horizontal. However, the slow line is occasionally punctuated by large avalanches that dip into 
the negative: curiously, the interface seems to tilt in the direction opposite that of rotation every once in a while. 
Figure 03 looks at one of these dips more closely. Apparently, even with the constant small avalanches that we see in 
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FIG. 10: Still frames depicting the momentum flow in our pile in the initial stages of rotation, in gravity's frame of reference. 
The darker arrows represent sites in the pile, while the lighter arrows are the very low-density sites above the pile. The length 
of each arrow is proportional to the momentum, pv, of that particular site. For legibility, each arrow is actually the average of 
four lattice sites. 



the bulk angle in the slowly rotating case, there is still a build-up of potential energy that must be released by these 
larger events. 

All of our tests so far show a nonzero dynamic angle of repose, but to rule out the possibility that these angles arc 
due to viscosity or inertia, it is best to find how the average angle (bulk or surface) depends on the period of rotation 
T, and from this relationship extrapolate to T — > oo. With this in mind we measure the average bulk and surface 
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FIG. 11: The number of close-packed sites for a cylindrical system rotated at a rate of T = 200. Plotted against this (and on a 
different t/-axis) is the bulk angle of the system, for comparison. Notice that the bulk angle does not stop its ascent until after 
the number of close-packed sites have begun to decline, and the pile has dilated. 




FIG. 12: The bulk angle over three revolutions for two different speeds of rotation. The horizontal axis is again scaled to show 
the angle through which the container has rotated, rather than the time taken. 



angles over time for several rotation speeds, and, in Figure ^ plot them versus the rotation speed 1/T. 

The first thing to note in these plots is that neither angle is going to zero in the limit of small angular velocities: 
there is a definite non-zero angle of repose in our system. This angle, about 2.5°, is quite small compared to experiment 
where typically one finds angles of repose on the order of 30° |2^, 23 1. 

It is interesting to fit the angles in Figure |lj to a power law 

= eo + cT'^. (26) 

Such fits are shown in the figure, with the parameters given in Table |[ The intercepts are positive, confirming a 
non-zero angle of repose. The fit to the surface angle is not quantitative, but for the average bulk angle the power 
law fits well, with an exponent of —1.4. 

These data points come from single runs, and averaging over an ensemble of initial conditions may make the fits 



TABLE 1: We fit the angles of repose in figure ^ to the model 6 = Oo + cT" . 

c a 00 



bulk angle 10700±4100 -1.35±0.07 2.49±0.15 0.148 
surface angle 83.7±105.0 -0.555±0.275 2.51±1.12 0.780 
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FIG. 13: An instance in the slow {T — 1600) rotation of the pile where the bulk angle dips below zero. The four flow diagrams 
show the momentum of the pile ai t — 1320, 1340, 1360, and 1380. The mechanism by which the angle drops is a large vortex 
in the pile which give individual grains a large momentum to the left. Notice how calm the pile is before and after this event: 
certainly this is intermittent behavior. 



10 
9 
8 
7 
6 
5 
4 
3 
2 
1 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
100/T 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 
100/T 



FIG. 14: The average surface (left) and bulk angles are the center points in each column of their respective plots; the outer 
points delineate one standard deviation above and below. Each data point represents one run of the system, consisting of the 
last two of three complete revolutions; the first turn was thrown out to diminish initial effects. The lines are power law fits. 



more quantitative. 

Experimentally, Rajchenbach finds that the angular velocity fl goes as 

n^id- (27) 

where m — 0.5 ±0.1. This corresponds to 

- T-2±0-4 (28) 

Of our two measures, the bulk angle comes closer to matching experiment, but there is still some discrepancy that 
needs to be addressed. 

Finally, to compare the surface and bulk angles with each other. For high-speed rotations, the bulk angle is higher 
on average than the surface angle, reflecting the plume of material that creeps up the side of the container (which 
is accounted for in a bulk calculation, but specifically excluded from the surface angle). The surface angle actually 
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seems to level off at high speeds. At lower speeds, the bulk angle actually dips below the surface angle, suggesting 
that the plume has disappeared (as is seen in Figure ^3|). Notice that the fluctuations in the surface angle are larger 
than those in the bulk angle (the former depending on fewer lattice sites and thus more volatile), and that in both 
angles these fluctuations get smaller as one reduces the speed of rotation (but do not seem to go to zero). 



VIII. DISCUSSION 



We present evidence that one can create a nonlinear hydrodynamical model for granular materials that depends 
only on density and momentum fields. Qualitatively, our model demonstrates many key features of sand, including 
a sharp interface, a relatively uniform density, a nonzero angle of repose, and a metastable structure. The method 
allows us to follow the pile from creation forward, in static and dynamic situations, and the model generally behaves 
like a sand pile. 

There are several ways in which the model could be improved. The angles of repose seen here are too small when 
compared to experiment, and the way that the angles scale with rotation speed is at odds with Rajchenbach's findings. 
This may not be a problem since we have not yet looked at the variation of the angle of repose with the parameters 
characterizing the model. An explanation for the difference may also lie in the fact that our rotation probe differs 
from the typical experimental method of rotating a sandpile. In our simulations, where we rotate gravity, every 
particle feels the external force directly and immediately. In experiment, where the container is rotated, the external 
force must be transmitted inward from the boundaries. This difference may be enough to account for the discrepancy 
between simulational and experimental outcomes. One may be able to mimic the rotation of the container in our 
model by introducing centrifugal forces into the system. Such forces would have a magnitude of 

Fc = pr i^ — 

where r is the distance from the center of the container. Our round container has a diameter of 100 units, so r < 50, 
while p 1 and the fastest rotation speed we use is T = 200; thus these centrifugal forces would have a magnitude of 
0.05, which is small (though not negligible) compared to the main acting force of gravity, which is of magnitude roughly 
equal to 1. Thus centrifugal forces may provide some quantitative effect, but in our initial, qualitative presentation 
here we deemed it unnecessary to include them. 

Another way of improving the model is to remove the constraint that the loose-packed regions have a single fixed 
density. There are many non-optimal, metastable ways to pack particles together, having a range of different densities. 
One solution may be to allow the position of the loose-packed well to fluctuate slightly at random through the pile. 
There are several plausible ways of implementing this idea, which we intend to pursue. 

We have only begun to extract information using our model; there are other probes we can use to perturb our system. 
Shaking can cause the pile to slowly settle into a denser state; we hope to flnd the logarithmic time dependence seen 
in compactification experiments pO| . Applying pressure to the system may allow us to investigate force propagation 
in the pile, and determine the nature of stress chains. It should be possible to modify the model to depict a pile 
made up of two or more types of particles, to investigate the phenomena of unmixing and the Brazil nut effect in a 
hydrodynamical setting. The flexibility of the fluctuating nonlinear hydrodynamical approach gives us a wide range 
of avenues to investigate. 

Appendix: Specification of the Free Energy Density 

The flrst potential in Figure |^ is made up of four terms: 

fa{p) = /clump(p) + /largo(p) + /nogativc(p) + /ontropy(/5) (29) 

where 

/clump = ^2^oP^ (30) 

/large = Be^^P"~^^ (31) 
/negative ^ UiC . (32) 

The flrst term models the clumping behavior of sand mentioned in property |l| in section II. The second term models 
the ultimate incompressibility of the sand grains: B is chosen so that the minimum of the well is at p = 1. The third 
term prevents the densities from becoming negative by putting a barrier at zero density. 
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We have not given the definition for /entropy yet. Our original intent was to have this term model the behavior of 
the very dilute gas that exists above our sand pile. So that the low-density regions show a Boltzmann distribution, 
we used the standard gas entropy term 

/entropy = ^iP / Po) - p) (33) 

with Po = 0.05. However, this was a source of numerical instability, and since our focus was the pile and the 
high-density regions, we replaced this with a simpler term, 

/entropy = Ap. (34) 

The reason we did not eliminate the term entirely was that the original term created a shallow minimum around 
p — 0.05, and for consistency we decided to keep that minimum there which the linear term allows us to do. 

Clearly, this model has a lot of parameters, and one part of our future work will be to find optimal values for 
these parameters. Our current choices were selected because they gave realistic results: we set uq = 2, A = 40, 
Ui = 0.6, Ai = 400, and A = 0.2. To satisfy the requirement that the minimum of the potential be at p = 1, we set 

The barrier introduced in Figure plb is described by 

/barrier ^Wbe-^'^t"-"')', (35) 

where Ub = 1.5, pb = 0.99, and k = 10^. The two wells in the final potential are described by 

Uon,,=-Une-i''^P-P"^\ (36) 
where pi = 0.98 and p2 — 1.01, ui = 0.25 and U2 = 0.4, and k is the same as in the expression for the barrier. 
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